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INTRODUCTION 


The  periodic  aerodynamic  loading  of  helicopter  rotors  leads  to  undesirable 
vibration  levels  throughout  the  aircraft.  These  vibrations  result  in  pilot  fatigue, 
passenger  discomfort  and  greater  maintenance  and  operating  costs.  For  military 
helicopters,  vibration  also  causes  difficulty  in  aiming  sights  and  pointing  weapons. 
A  reduction  in  operating  vibration  levels  would  therefore  be  of  benefit  to  the 
military  and  commercial  helicopter  community. 

The  helicopter  vibration  problem  can  be  attacked  by  two  alternate  strategies. 
In  the  first,  the  helicopter  or  components  of  it  are  isolated  from  the  rotor  so  as  to 
reduce  the  forces  or  the  vibrations  transmitted.  The  aerodynamic  loading  remains 
unchanged  by  such  response  attenuation  devices.  The  second  strategy,  source 
alleviation,  attempts  to  remove  the  periodic  airloading  through  an  azimuthal 
alteration  of  the  blades'  lifts.  Response  attenuation  devices  have  been  investigated 
extensively  but  are  largely  unsatisfactory  due  to  their  weight  and  insufficient 
robustness.  Source  alleviation  methods  have  been  employed  and  examined  in  a 
greater  number  of  studies  and  have,  in  general,  yielded  promising  results.  Many  of 
these  approaches  have  included  on-line  system  identification  so  that  optimal 
vibration  reduction  can  be  achieved  as  the  helicopter  plant  changes  \vitU  Right 
condition  [  1  — 5S] .  In  this  report,  an  adaptive  source  alleviation  controller  is 
formulated  m  terms  of  a  "ew  input— output,  model.  For  any  detail  on  this  research 
effort,  please  see  F  fence  tiT 


CHAPTER  ONE:  THE  IMPULSE  B ESPOXSE  METHOD  AM)  ESTIMATION 


1.1  Introduction 

The  impulse  response  method  seeks  to  find  a  time  sequence  of  discrete 
control  inputs  to  the  blades  that  will  result  in  highly  reduced  fuselage  vibrations. 
The  method  relies  upon  the  quasi— periodic  nature  of  the  uncontrolled  fuselage 
vibrations  and  the  plant  dynamics.  This  yields  a  predictable  plant  where  previous 
measurements  can  be  used  to  identify  the  periodic  time  sequences  describing  the 
uncontrolled  vibration  and  the  impulse  response.  The  impulse  response 
characterizes  the  plant  dynamics.  If  the  plant  is  linear  (which  it  should  be  for  small 
control  amplitudes)  then  the  control  induced  output  is  just  the  convolution  of  the 
impulse  response  and  the  control.  Discretized,  this,  as  shall  bo  shown,  takes  the 
form  of  matrix  relationship  relating  the  control  time  sequence  and  the  vibration 
lime  sequence.  (If  the  flight  condition  is  changing  slowly,  the  plant  is  periodic,  and 
the  matrix  has  repeating  blocks  whose  length  is  one  period  of  the  plant  along  its 
diagonal.)  The  matrix  is  identified  by  an  estimator  along  with  the  uncontrolled 
vibration  in  some  formulations.  The  estimates  are  then  used  by  an  optimal  control 
law  to  compute  a  control  sequence  for  the  next  plant  period.  For  this  research,  the 
vibration  reduction  was  to  be  obtained  only  on  the  vertical  axis.  This  allows  the 
same  control  sequence  to  be  fed  to  each  blade.  (For  multi-axis  vibration  reduction 
the  blades  must  execute  different  control  sequences  simultaneously  -  see  Reference 
63.)  Thus,  at  any  time,  all  blades  are  executing  the  same  deflection.  The  period  of 
the  plant  in  this  case  is  one  quarter  revolution  which  in  this  study  is  called  one 
cycle.  Therefore,  all  control  input  sequences  are  harmonics  of  10.  The  controls 
commanded  by  this  system  could  be  carried  out  through  a  swashplate  although  this 
actuator  would  require  a  higher  actuator  bandwidth  than  on  present  helicopters. 
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For  the  multi— axis  case,  however,  individual  blade  control  must  be  used.  In  either 
case,  the  mechanics  of  blade  control  are  not  essential  to  the  derivation  of  the 
vibration  control  system  and  are  not  considered  here. 

1.2  Impulse  Response  Method 
For  anv  linear  system, 

z(t)  =  A(t)z(t)  +  B(t)f(t) 

y(t)  =  C(t)z(t)  +  D(t)f(t)  +  v(t) 

the  response  y(t)  to  a  given  forcing  function  is 

.v(t)  =  C(t)  (*(t,t0)z(t0)  +  /  4>(t.,r)B( r)f( r)dr]  +  D(t)f(t)  +  v(t) 

Co 

and  for  a  system  where  the  initial  conditions  have  been  damped  out  and  the  forcing 
term  is  the  combination  of  a  disturbance  d(t)  and  a  second  derivative  control  u(t) 
the  response  takes  the  form 

t, 

y(t)  =  D ( t ) d ( t )  /  C(t)4>(t,r)B(r)d(r)dr  +  I).,(t.)u(t) 

lo 

t 

+  /  C(t.)4>(t  ,r)B.,(  r)u(  r)dr  +  v(t) 
t 

o 
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If  the  disturbances  d(t)  and  the  impulse  response  4>(t,r)B(t)  were 
predictable,  then  an  open  loop  control  could  be  found  to  minimize  a  performance 
index  of  y  and  u.  The  response  can  be  discretized  by  replacing  the  integrals  with 
summations 

y(tj)  =  Dltjjdftj)  +  £  C(tj)*(tjltk)B(tk)d(tk)  +  D2(tj)u(tj) 

+  £  C(t.)$(ti,tk)B.,(ti)u(ti)  +  v(t.) 
k=o  J  J  K  -  J  J  J 

which  can  be  rewritten  as 


x<y = xo<y  +  w(tj,tk)u(‘k) 

y(V  -  x(v +  v(tj> 

or  in  matrix  form. 


If  the  disturbance  responses  x0(tj)  and  the  control  impulse  response  w(tj,tj.) 
(i.e.  the  weighting  function)  are  predictable  then  an  open  loop  control  could  be 
found  to  minimize  a  performance  index  of  y  and  u. 


1.3  Application  to  Helicopter  Vibration  Reduction 

For  a  helicopter  in  a  steady  flight  condition,  the  induced  vibration  can  be 
considered  periodic 


xo^in+j)  xo^in+n+j^ 

where  n  measurements  are  made  per  cycle,  which  is  one  quarter  of  a  revolution  for  a 
four-bladed  rotor  when  only  a  single  vertical  accelerometer  is  assumed.  Thus,  if  the 
induced  vibration  history  can  be  determined  for  one  cycle,  it  can  be  predicted  to  be 
the  same  for  the  next  few  cycles.  If  the  effect  of  small  blade-lift-controls  are 
negligible  alter  p  cycles  then  the  matrix  equation  becomes 


+ 

w 

1 

■  II 

u<lin+l> 

zs~ 

II 

xo^in+l^ 

x(t.  ,  ) 
v  111  +  11' 

+ 

xo^in+n^ 

x. 

1 


l—i 


+  V 
1 


(1.1) 


where  Rj_]  is  an  impulse  response  matrix  relating  cycle  i-1  input  to  cycle  i  output. 

If  the  control  is  nearh  constant  from  one  cycle  to  the  next,  so  that 
u(tjn+j)  ~  u(tjn+n+j),  then  the  matrix  equation  can  be  written  for  cycle  i 
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x.  =  Tu.  +  x 

1  11  Oj 

(1.2) 

y-,  =  xi  +  vi 


T.  =  E  R.  . 
1  1=0 


where  T-  is  the  "wrapped  around"  impulse  response  matrix.  Throughout  the  rest  of 
this  report  this  shall  be  cailed  simply  the  impulse  response  matrix. 

The  impulse  response  method  thus  seeks  to  find  vector  m  such  that  a 
performance  index  of  the  vibration  vector  x^  is  minimized.  Note  that  a  more  general 
regulator  algorithm  can  be  obtained  by  using  Eq.  (1)  without  the  nearly  constant 
control  assumption.  Such  an  algorithm  would  be  more  complex,  but  would  penult 
treatment  of  transients  generated  by  the  control  inputs  more  effectively. 

This  development  may  be  generalized  such  that  u  is  a  vector  of  m  discrete 
blade  inputs  evenly  distributed  throughout  the  blade  azimuth,  y  is  ameasure  of  the 
vector  x  of  u  discrete  values  of  an  observable  response  (typically  accelerations) 
corrupted  by  noise  \\  T  is  a  nxm  transfer  matrix  composed  of  impulse  responses, 
and  xQ  is  a  n— vect  >r  of  discrete  values  of  the  uncontrolled  vibrations.  By  the 
Sampling  Theorem,  the  number  of  samples  per  cycle  taken  of  the  vibration  (n)  and 
used  to  reconstruct  the  control  signal  (m)  must  be  greater  than  twice  the  highest 
harmonic  frequency  of  vibration  to  be  reduced.  This  rate,  however,  should  be  at 
least  four  times  Nyquist  for  proper  performance.  In  this  study,  the  highest 
frequency  in  the  uncontrolled  vibration  was  Stt  yielding  a  Nyquist  frequency  of  16 
samples/  revolution,  or  4  samples/cycle  (one  cycle  equals  one  quarter  revolution). 
The  rate  used  here  was  S  samples  per  cycle;  therefore  m  =  n  =  8  (twice  Nyquist). 
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In  the  following  development  the  subscript  i  represents  the  current  cycle  count; 
however,  in  the  interest  of  simplicity,  the  indices  are  only  used  on  expressions  which 
represent  relationships  between  different  cycles. 

Eq.  (1.2)  can  be  rewritten  in  a  more  convenient  form 

x  =  On  (1.3) 

where 

«=lx0|T] 

and 

r,  i  T,T 
u=  [1  |  u  ] 

Eq.  (1.2)  is  in  global  model  form.  It  can  be  written  in  local  model  form 
through  simple  subtraction 


where  it  is  assumed  that  Tj  j 
into  the  same  form  as  ( 1.3) 


T:  and  x  =  x  .  Equation  (1.4)  can  be  put 
°i+ 1  l 


Ax=  On 
0=\ T] 
u  =  Aii 


s 


The  local  model  requires  less  estimation  to  be  performed  for  the  optimal  control 
law.  It  is  also  felt  by  HHC  researchers  to  better  represent  the  globally  nonlinear 
relationship  between  blade  input  and  fuselage  vibration.  For  small  changes  in 
con'rol  from  a  reference  blade  deflection,  the  change  in  control  force  (and  hence 
counter  vibration)  produced  will  be  linearly  related.  However,  since  the 
effectiveness  of  changes  in  blade  deflection  varies  with  blade  deflection,  this  linearity 
is  only  valid  locally  (that  is,  close  to  the  reference  deflection).  Equation  (1.4) 
expresses  only  a  locally  linear  relationship,  not  a  globally  linear  one  as  (1.2). 

1.4  Circulant  Structure  of  the  Impulse  Response  Matrix 

If  the  plant  was  not  a  time  varying  system,  each  column  of  the  square 
impulse  response  matrix  would  be  an  identical  permutation  of  the  first  column: 


To  Tj 
T3  T2 


'u-1  ’ 
Tn  • 
T1 


..  Tr 


Tn  Tn— 1  Tn-2 


Thus,  the  reaction  of  an  impulse  in  control  at  any  point  during  a  cycle  would  result 
in  the  same  response.  A  matrix  of  this  form  is  called  a  circulant  matrix  and  has 
many  special  properties.  Of  particular  interest  is  that  every  circulant  matrix  is 

x 

diagonalizable  by  the  Fast  Fourier  (F)  and  inverse  Fast  Fourier  transform  (F  ); 


therefore  if  T^  =  circ  (Tj.  T.,,  ...  l’n) 


T  =  F  AF 


x  —I 
F  =  1-'  1 

A  =  diag  (A, . A|() 
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and 

Fx  =  F(F  AF)u  +  Fxq 

f>j  .  r-j  rv i 

x  =  A  u  +  x0 

where  the  tilde  indicates  the  FFT  of  the  time  sequence  vector.  Thus,  for  the 
stationary  plant,  the  circulant  impulse  response  matrix  expresses  the  spectral 
response  of  the  transfer  function  of  the  plant  through  its  eigenvalues. 

If  the  vertical  helicopter  impulse  response  at  each  blade  azimuth  can  be 
idealized  as  the  same  except  for  an  azimuth  dependent  magnification  factor  each 
column  of  the  impulse  response  matrix  will  be  a  permutation  of  the  first  column 
multiplied  by  a  scalar.  The  magnification  factors  are  then  just  an  expression  of  the 
changing  effectiveness  of  the  blade  (i.e.,  coefficient  of  lift)  throughout  the  azimuth. 
(The  force  produced  on  a  blade  at  any  point  in  the  azimuth  by  an  incremental 
change  in  blade  pitch  is  a  function  of  the  azimuth  dependent  blade  pitch.)  This 
idealization  ignores  all  other  azimuth  dependent  dynamics  (such  as  azimuth 
dependent  blade  damping).  The  author  has  not  seen  any  research  that  suggests  that 
the  effects  of  the  ignored  dynamics  will  be  substantial.  Therefore,  the  helicopter 
impulse  response  matrix  can  be  considered  decomposable  into  a  circulant  matrix,  C, 
and  a  diagonal  matrix.  B: 

T  =  C  B 

T 

C  =  circ(cr  c2.  ...  cn) 

B  =  diag  (bj.  b.,,  ...  b  ^ 

Note  that  this  decomposition  is  not  necessary  for  the  application  of  the  impulse 
response  method  to  the  helicopter  vibration  problem.  Note  also  that  this 
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decomposition  requires  that  the  control  be  reconstructed  from  the  same  number  of 
samples  as  are  measured  per  cycle  (i.e.,  in  =  n).  This  presents  no  problem  since  the 
control  should  have  a  bandwidth  as  large  as  that  of  the  uncontrolled  vibration. 

Applying  the  circulant  diagonalization  property  to  the  global  model  yields 

x  =  C  B  u  I  xy  (1.5) 

Fx  =  F(F*AF)(Bu)  +  Fxq 

x  =  A(flu)  +  xQ 

Thus,  each  vibration  harmonic  is  the  sum  of  the  corresponding  uncontrolled 
vibration  harmonic  and  the  corresponding  forcing  (i.e.  Bu)  harmonic  attenuated  and 
phase  shifted  by  a  complex  eigenvalue  of  the  circulant  matrix.  This  comes  close  to 
relating  the  impulse  response  method  to  the  HHC  method;  HHC  relates  the 
harmonics  of  vibration  (x)  to  harmonics  of  control  (u)  via  a  full  matrix.  The 
impulse  response  matrix  formulation  can  be  put  into  this  form: 

A(F  B  F*)  Fu  +  xQ 
A  Hc  u  +  xQ 

II  =  FBF*  =  circ  (h. .  h0,  ...  h  ) 
c  M2  n ' 

where  11^  is  a  Hermetic  circulant  matrix. 


1.5  Computer  Models  of  the  Plant 


The  helicopter  model  employed  for  the  computer  simulations  to  be  presented 
is  based  upon  equation  (1.5) 

x=CBulxQ 
y  =  x  +  v 

T 

C  =  circ  (cp  c2,  ...  cn) 

B  =  diag  (bp  b.,,  ...  bn) 

where  x,  x  .  y,  u  and  v  are  8-vectors,  and  C  and  B  are  SxS  matrices.  The  8 
elements  of  xQ  are  determined  by  4  parameters:  a  4 Cl  amplitude,  a  4f 1  phase,  an  Sfl 
amplitude,  and  an  8 ft  phase.  The  8  non— zero  elements  of  B  are  determined  by  3 
parameters:  a  bias,  a  sine  function  amplitude,  and  a  sine  phase.  This  models  the 
one  peak  per  cycle  combined  blade  lift  sensitivity.  The  elements  of  the  circulant  C 
matrix  are  determined  by  6  parameters  (these  determine  the  first  column  of  the 
matrix,  the  wrapped— around  impulse  response),  3  for  each  clamped  vibration  mode 
in  the  impulse  response:  an  amplitude,  a  frequency,  and  a  damping.  The  impulse 
response  is  computed  for  3  cycles  and  overlapped  to  produce  the  first  column  of  the 
circulant  matrix.  Thus,  13  parameters  determine  all  72  elements  of  the  plant  model 

<v  T>- 

1.6  Kalman  Fiber  Akorit  Inn 

To  adapt  to  the  changes  in  the  plant,  the  regulator  must  incorporate  an 
identification  algorithm.  The  filter  developed  here  parallels  that  used  by  Molusis  et 
al.  [34],  and  Bryson  and  Ho  [59].  The  0  matrix  varies  from  cycle  to  cycle  as  the 
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result  of  disturbances  in  the  atmosphere,  flight  control  input  from  the  pilot,  or  self 
induced  disturbances,  such  as  interactions  between  the  blade  wake  and  the  tail 
rotor.  These  variations  are  modeled  by  the  random  walk. 


"i+i  =  °\  +  w«.i 


where  j  is  a  matrix  of  random  numbers 


w0  =  i  wx  1  WT 


w,,  =  w 


global 


9  -  l  "T 


Here  w  is  a  vector  of  n  random  variations  in  the  basic  vibration  vector  xq 
with  standard  deviation  a  and  is  an  nxm  matrix  of  random  variations  in  the 
transfer  matrix  T  with  standard  deviation 

The  value  of  0  differs  from  the  minimum  variance  a  po. staiori  estimate  0  by 
an  error  matrix  r/,, 


0- 0 


v  0  =  l  Vx  I  V T 

v0  =  [  VT  ] 


global 


in  which  rj^  is  a  vector  representing  the  error  in  the  uncontrolled  vibration  estimate 
vector  x  .  while  rjj  is  a  matrix  representing  the  error  in  the  impulse  response 
matrix  estimate  'I'.  The  minimum  variance  a  priori  estimate  of  +  j  is  #  +  1  and  is 
equal  to  the  previous  a  posteriori  estimate  0 .  so  that 
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Equations  (1.2),  (1.3)  and  (l.G)  can  now  be  combined  into  the  state  propagation 
equation  which  is  written  as 

Y=0H  +  (  (1.7) 

with 

Y  =  (y  i  D\ 

II  =  [  I!  !  1  1 

<  =  [y  I  <«] 

^.i+1  =  Ji+I  _  ®i+l 

=  lo. i  -  Ww 

1.7  Best  Linear  I'ubiased  Prediction 

At  the  end  of  the  it h  cycle,  several  quantities  must  be  updated  so  that  a  new 
control  vector  ^  can  be  found.  The  a  priori  estimate  #  +  j  is  equal  to  0 t  which  is 
obtained  from  the  fnea.surement  vector  of  tin1  helicopter  dynamics  in  response  to 
the  input  of  u-  A  linear  relationship  is  assumed 

0  -  Y  M 

The  requirement  of  an  unbiased  estimate.  Y.(0)  =  0.  yields  the  condition 

II  M  =  I 


u  —  [1  |  u1]1  global 
=  Au  local 


i 


1-1 


The  best  estimate  is  obtained  by  minimizing  the  error  covariance  matrix,  Pj. 


subject  to  this  condition.  The  matrix  P.  is  defined 


T 


pj  =  E(,0.’(d 

J  J  J 


where  E{-}  is  the  expected  value  operator  and  77^  is  the  j'th  row  of  the  estimate 

j 

error  matrix.  This  minimization  yields  the  Kalman  filter  algorithm  for  the  estimate 
of  j'th  row  of  0 


P.  =  P.  -  P.  u.KT  +  W. 
Jj+1  Jj  Jj  1  J  j  J 


(l.S) 


D, 


T 


•‘i+l 


0-  =  0  +  (y-  -  0 ■  u.)K  . 

1.  1.  VJt.  ,-r 

j,  Jj  j,  j,  j  , 


Here.  K.  is  the  Kalman  gain  vector 

1  c> 


K.  .=  P.  u./(uf  P  ti  +  V. 
1.  1.  — 1 '  v“  1  1—1 


Jj  Jj-P'-i  Jp  J 


while 


=  K{  v  v}  =  a 

.1  .1  y  v 


global 

local 


W,  =  E  10.  -0.  )'{0. 

1  h-t  1  Ji  Ji  t  1  ■  i 


‘)  ) 


0:  )}  =  diagjrT.  a.j. . rr  j  j  global 


=  (iiag  { rrT- .  <T- . 


e.,. ) 


local 


0 


and  Pj  is  the  a  priori  estimate  error  covariance  matrix 


pj  = E  -  °i '  q  -  "j'l 


Therefore,  for  the  estimation  of  each  row  of  0 ,  an  nxn  matrix.  Pj,  must  be  updated 

(and  stored)  and  an  nxl  vector,  K-,  must  be  computed.  These  are  dependent  on  the 

noise  covariance  on  the  j'th  measurement,  Vh,  and  the  process  noise  covariance 

matrix  for  the  j'th  row  of  0 ,  W-.  There  is  no  reason  to  assume  in  advance  that  these 

J 

quantities  will  be  different  for  different  rows;  and.  considerable  computational 
advantage  can  be  achieved  by  assuming  them  equivalent.  In  this  case,  the 
individual  row  equations  (i.S)  can  be  collated  into  a  single  matrix  update  equation 
for  the  estimate  and  a  single  equation  for  the  covariance  matrix 


Ji+i = +  <>i  -  j»p  kT  <l9> 

Pi+l=Vb!li  KI  +  W 

K.  =  F.  u./(uT  F.  u.  +  V) 

i  l  -y  v- 1  l  -i  ' 

\7  =  a~  global 

=  2  a“  local 

v 

W  =  diag  { a”,  . a.p}  global 

=  diag  { i7 -p . <7,‘ j.}  local 
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> 


I 


> 


I 


I 


! 


I 


* 


This  reduces  the  computational  burden  n— fold.  While  in  a  typical  flight,  the  actual 
variation  in  each  row  of  <7  will  be  different,  it  is  difficult  to  prescribe  the  row  process 
noise  covariance  matrices  in  advance.  Therefore,  this  assumption  is  also  practical. 
This  will  result  in  a  degradation  in  the  estimation  performance  in  comparison  to  the 
optimally  tuned  individual  row  filters. 

The  Kalman  filter  algorithm  (1.9)  recursively  updates  a  measure  of  its 

uncertainty,  P,  and  the  a  priori  estimate,  0.  The  change  in  the  estimate  of  i'th 

column  produced  by  the  update  is  the  residual  vector  (yj  —  #u.)  multiplied  by  the 

j'th  element  of  die  Kalman  gain  vector  K^.  The  residual  vector  is  the  only 

expression  of  the  error  in  the  estimate  available  to  the  filter;  it  is  the  measured 

output  vector  less  the  expected  output  vector  based  upon  the  estimate.  The 

Kalman  gain  vector  expresses  the  usefulness  of  this  information  to  the  updating  of 

each  column  of  the  estimate  matrix.  The  larger  the  magnitudes  of  the  elements  of 

the  Kalman  gain  vector  the  greater  the  change  in  the  estimate  will  be  and  the 

—  T 

greater  the  reduction  in  the  covariance  ellipsoid  via  the  PjiijK  ^  term.  Quite 
simply,  if  the  residual  information  is  very  useful,  then  the  expected  error  of  a  new 
estimate  based  on  it  should  decrease  from  the  previous  estimate.  (Note  that  the 
usefulness  of  the  residual  is  inversely  proportional  to  the  expected  measurement 
noise  covariance.  V.)  Of  course  any  change  in  the  the  value  after  the  measurement 
is  taken  will  increase  the  estimate's  inaccuracy,  hence  the  addition  of  the  process 
noise  term.  W.  to  the  previous  covariance. 

If  the  plant  is  stationary  and  the  measurement  noise  is  small,  the  Kalman 
filter  estimates  will  converge  to  the  true  values  in  approximately  m+1,  for  the 
global  model,  or  m.  for  the  local  model.  steps.  Also,  if  the  plant  is  stationary  the 
Kalman  filter  algorithm  is  stable  if  the  covariance  remain  bounded.  The  greatest 
problem  in  employing  the  Kalman  filter  is  the  extensive  tuning  (that  is  adjustment 


* 


17 


of  the  measurement  noise  and  process  noise  covariances)  required  for  near— optimal 

performance.  In  this  study,  it  is  assumed  that  the  measurement  noise  standard 

deviation  will  be  known  very  accurately;  therefore,  tuning  refers  here  to  the 

adjustment  of  the  process  noise  terms,  which  are  the  filters  representation  of  the 
2  ,  2 

true  covariances,  a  and  <7™. 

’  x  i 


1.8  Kalman  Filter  Covariance  Stability 

The  Kalman  filter  equations  (1.9)  can  become  unstable  via  growth  in  the 
covariance  ellipsoid  if  the  control  fails  to  yield  information  along  axes  of  the 
parameter  space.  For  example,  if  all  the  elements  of  the  control  vector  are  always 
zero  except  the  k'th  which  is  always  one,  then 


pk,  pk. 
=  P,  +  W  +  — ! - 


1 


fh.  +  V 
kkj 


(1-10) 


where  P^  is  the  k'th  column  of  P  and  is  the  k'th  element  of  Pj..  Then, 
equations  for  the  elements  of  P  can  be  written: 


P 

P, 

P 


11.,. 

=  pn. 

+  wn 

1+1 

l 

kk 

=  Pkk 

+  wkk 

l  +  l 

l 

lki+l 

=  pik. 

i 

-p>", 

-plk./<pkk.  +  v> 

i  i 

-  pkk./(pki, +  v> 
i  i 

pkk./<pkk.  +  v» 


(l.m 


The  second  equation  converges  to  the  stable  root 


pkk  =  I  (wkk  +  ^  wkk +  4Wkkv  i  > 0 

oo 


As  this  occurs,  the  last  equation  uniformly  converges  as 


p 

M<k 

Pi,.  =P,|-  (l-^-25 - 

1  i+I  iki  P,.  +  V 


) 


kk 


oo 


to  zero  (P^  =  0)  since  0  <  (1  -  (P^  +  V)  *)  <  1.  Therefore,  the  first 

00  00  voo 

equation  becomes 


n.  =pn.  +  wii 

i+i  iXi 


ik 


Apkk  +v) 


00 


00 


or 


+  W 
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(1.12) 


as  the  second  and  third  equations  approach  their  roots.  It  is  easy  to  see  from  (1.12) 
that  this  covariance  will  go  unstable  under  this  control.  Furthermore,  it  can  be 
shown  in  the  same  manner  that  all  Pjj.  j  #  k,  will  go  unstable  under  this  control. 
The  covariance  ellipsoid  is  growing  along  all  axes  except  the  k'th  along  which  the 
filter  is  receiving  information.  This  may  seem  like  an  unusual  case;  after  all,  the 
control  is  not  likely  to  be  zero  in  every  element  but  one.  However,  if  the  control  is  a 
constant  vector,  then  through  an  orthogonal  transform  L  we  can  rotate  the  axis 
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i 


UL  =  Lu  = 


0  0  ...  0  1  0  ...  0] 


(In  general,  the  magnitude  of  t.he  constant  control  can  always  be  one  if  the 
measurement  and  process  noises  are  also  scaled.)  Then,  by  applying  the  transform 
to  the  covariance  equation, 


LPj+iLT 


Lf.LTtLWLT-L{fr)[AT)LP  LT 

1  (u1!,1  )(L  P  L1 )  (Lu)  +  V 


(which  is  geometrically  equivalent  to  rotating  the  covariance  ellipsoid  in  the 
n— space)  and  the  substitution 


fl  =  lplt 


Lu 


the  recursive  equation  for  P  is  converted  to  the  form  of  (1.10). 

An  examination  of  the  case  where  one  element  of  the  control  is  zero  (and  the 
others  may  vary)  is  useful.  The  covariance  update  equation  (1.9)  can  be  partitioned 
as 


Pa 

1 

-O  1 

ICu  i 

'  Pa 

pb ' 

’  Wa 

0 

~*r 

P 

nn 

— 

-T 

i>  1 

1  b 

p 

nil 

+ 

0 

w 

i  +  1 


pa 

pl> 

ua 

u l  |  0 

'  ! 

P 

a 

pb 

Pb 

p[in  . 

i 

0 
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i 

: 
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+  V 
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to  yield 


P  u  uT  P 

_  a-  a-  a-  a- 

Pa  =  Pa  +  Wa  - 

i+l  ‘  pfl  +  v 

a  i  Ai  ai 


b- 


i+l 


^Pa.ua.ua.  > 

T _ I _ 1 _ I _ 

uIPa.V+  V 

1  i  i 


b; 


(“l.  Pb-)2 


nn. 


i+l 


P  +W  ----- - - J - 

nn.  nn  T  5  ,  ,, 

1  Pa  Ua  +  V 

a  i  ai  ai 


The  bracketed  matrix  in  the  second  equation  of  (1.13)  has  eigenvalues 


uaiPajuai 

Al  =  1-3^ 


>  0 


VW v 


Aj  =  1 


J  #'1 


Since  all  the  eigenvalues  are  less  than  or  equal  to  one,  any  norm  of  P, 

i+l 

increase  from  i  to  i+l  Therefore,  over  time  the  magnitude  of  the  vector 
decrease;  as  this  occurs,  the  last  equation  will  increasingly  resemble 


5  =  P  +  W 

nil.  ,  ,  nn.  nn 

i+l  i 


(1.13) 


cannot 
Pb  will 
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Thus,  with  one  element  of  the  control  always  zero,  one  axis  of  the  covariance 
ellipsoid  will  tend  to  grow  unbounded.  This  is  geometrically  equivalent  to  the  case 
where  the  control  is  constrained  to  lie  in  a  space  orthogonal  to  a  fixed  vector;  that 
is,  through  a  coordinate  transformation  one  axis  of  the  basis  can  be  aligned  with  the 
fixed  vector  yielding  a  representation  of  the  control  with  the  corresponding  element 
fixed  equal  to  zero.  Thus,  if  the  control  is  constrained  to  have  zero  mean,  then  one 
eigenvector  of  the  covariance  ellipsoid  will  grow  unbounded.  This  growth  in 
covariance  represents  the  increasing  uncertainty  in  the  filter  over  the  estimate  of  T. 
If  every  element  of  T  were  to  increase  by  the  same  constant,  then  this  change  would 
be  unrepresented  in  the  residual  when  the  control  is  constrained  to  have  zero  mean. 
While  this  behavior  is  exceedingly  unlikely,  the  covariance  equations  instability 
represents  the  chance  of  it  occurring  given  the  random  walk  description  of  the  plant 
model. 

These  arguments  can  be  generalized  to  say  that  the  number  of  bounded 
eigenvalues  as  the  recursion  approaches  infinity  is  equal  to  the  minimum  rank  of  the 
control  matrix 

U^=  u^,  u^+1  |  ...  u^+1  ...  ,  0  <  (  <  00  (1.14) 

over  finite  t  This  behavior  can  be  seen  in  Figures  1.1  through  1.3  which  show  the 
growth  in  the  eigenvalues  of  a  9x9  covariance  matrix  under  rank  9,  S,  and  7  control. 
(The  control  moves  mod  p  through  a  set  of  p  unit  vectors.)  Note  that  the  number  of 
'bounded'  eigenvalues  is  that  predicted  by  the  minimum  rank  of  (1.14). 

Since  in  any  sustained  flight  condition  the  optimal  control  may  be  almost 
constant,  it  is  likely  that  in  operation  an  on-board  Kalman  filter  will  often  produce 


uncertainties  much  larger  than  the  actual  estimate  errors.  This  would  especial1}' 
degrade  the  performance  of  cautious  regulators.  This  covariance  growth  can  be 
prevented  through  a  variety  of  means,  such  as  process  noise  turn-off,  random 
probing  control  and  covariance  trace  bounding. 

1.9  Weighted  Least  Squares  Algorithm 

The  moving  batch  least  squares  filter  algorithm  finds  the  estimate  that 
minimizes  the  sum  of  the  past  t  residuals  squared.  If  the  measurement  noise  was 
zero  and  the  plant  was  stationary,  this  method  would  produce  perfect  estimates 
with  any  batch  size  greater  than  or  equal  to  the  minimum  of  m+1  (global)  or  m 
(local).  If  the  plant  is  slowly  changing,  the  estimate  accuracy  is  increased  by 
weighting  the  most  recent  residual  squares  greater  (if  the  batch  size  is  not  equal  to 
the  minimum).  The  derivation  of  the  weighted  least  squares  (WLS)  algorithm  is 
straightforward.  Define  the  collections  of  i  measurement  vectors  and  control  vectors 
(where  the  subscript  B  stands  for  batch)  as 

YB  =  yi-l  |  yi-2  I  "•  I  yi-<  UB=  Vl  I  “i-2  I  -  I  «i-< 

Then  the  collection  of  residuals  is 


where  D-  is  a  prediction  of  #.  A  performance  function  that  is  the  weighted  sum  of 
the  f  residual  vectors. 
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J  = 


( 

v 

k=l 


^i- 


,  —  Du 

-k  1  l- 


-k)  ^k^i-k-^i-k) 


where  are  the  weights,  can  be  shown  equivalent  to  the  performance  function 


J  =  trace  (WB  [Yg  -  ^Ug] 1  [Yg  -  flUg]) 


WB  =  diag  (<pv  d>2,  ...,  <p£) 


Since  arguments  within  the  trace  parentheses  can  be  permuted  without  affecting  its 
value,  this  can  be  re— arranged  to 


J  =  trace  ([Yb  -  JjUJ  WB  |Yb  -  Sub]t) 


It  is  easy  to  show  that  when  the  first  variation  is  taken  of  J,  the  variations  in  the 
right  hand  side  may  be  taken  within  the  trace  parentheses  yielding 


dJ  =  trace  (dD{  {-Ug  Wg  Yg  +  Ug  Wg  Ug  flT} 


+  {  — Yb  Wb  U l  +  ^UbWuuT)  ^T) 


Setting  OJ  equal  to  zero  yields 


*iUBWBUB  =  ybwbub 
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Therefore,  the  weighted  least  squares  estimate  of  0 1  is 


^WLSj  “  YBWBUB  (UBWBUB) 


(1.15) 


This  same  result  may  alternately  be  shown  by  minimizing  the  n  performance 
functions 


Jk  =  <YB,  UB>  WB<YB.  ~\Ub) 

K  1  K  1 


T 


with  respect  to  0,  ,  where  the  k  subscript  on  Yg  and  #  refers  to  the  k'th  row  of  the 
i 

matrix.  This  yields  the  YVLS  estimate  of  0, 


^lJwLS  “  YB,  WBUB  ^UBWBUB^ 

1  K 


The  n  estimate  rows  can  then  be  collected  together  to  yield  (1.15). 

Note  that  the  control  batch  matrix  Ug  must  have  rank  m+ 1  for  the  global 
model  and  rank  m  for  the  local  if  equation  (1.15)  is  to  be  used.  Without 
constraints,  this  is  highly  likely  although  the  matrix  may  still  be  poorly  conditioned 
(especially  in  steady  flight  conditions  when  the  control  may  change  very  little). 
Care  must  be  taken  to  minimize  numerical  error  in  these  computations.  (A  special 
IMSL  routine  was  used  in  this  study's  simulations.)  Increasing  batch  size  will,  in 
most  cases,  result  in  a  better  conditioned  matrix  for  inversion.  For  a  stationary 
plant,  an  increase  in  batch  size  should  produce  greater  noise  suppression  and  better 
estimation.  If  the  plant  is  varying  quickly,  the  estimates  will  lag  the  true  values 


more  when  the  larger  batch  sizes  are  used.  Therefore,  generally,  a  batch  size 
slightly  larger  than  the  minimum  will  estimate  with  the  most  accuracy,  a  trade-off 
between  noise  suppression  and  filter  sensitivity. 


The  advantages  of  the  WLS  algorithm  are  that  it  is  always  stable  (since  it  is 
non— recursive)  and  it  needs  no  'tuning'  except  for  the  weighting  parameters.  For 
the  research  conducted  here,  the  weighting  scheme  employed  by  Jacklin  and 
Franklin  and  Powell  is  used  [58]: 


<K  =  (a) 


k— 1 


where  a  is  a  fixed  constant  between  zero  exclusive  and  one  inclusive.  One 
disadvantage  of  the  method  for  this  application,  is  that  a  zero  mean  control  batch 
can  never  be  of  rank  m+1  for  global  or  rank  m  for  local.  This  may  be  remedied  by 
adding  a  small  random  noise  to  the  control,  in  effect,  breaking  the  constraint 
deterministically,  but  not  stochastically. 

1.10  Circulant,  Form  and  Estimation 

If  the  true  plant  has  the  input/output  relationship  of  (1.5), 


x  =  G  B  u  +  x 
y  =  x  +  v 


T 

C  =  circ  (cr  c2.  ...  cn) 
B  =  diag  (bp  b.,,  ...  bj 


26 


then  only  3n  parameters  entirely  describe  the  plant  as  opposed  to  n“  +  n. 
Therefore,  for  a  stationary  plant  with  zero  measurement  noise,  the  unknown 
parameters  can  be  solved  for  with  only  3  measurement  vectors  taken,  as  opposed  to 
n+1.  It  is  easy  to  see  that  if  the  plant  is  described  in  this  form  the  estimates  of  the 
impulse  response  matrix  can  track  quick  changes  in  the  plant  (such  as  produced 
during  maneuvers)  with  much  greater  accuracy.  Unfortunately,  however,  since  the 
parameters  to  be  estimated  are  multiplied,  the  filtering  problem  is  now  nonlinear. 
Therefore,  either  a  linearization  needs  to  be  performed  so  that  an  extended  Kalman 
filter  (EKF)  can  be  used,  or  a  full  scale  nonlinear  filter  must  be  derived.  This 
second  course  of  action  is  exceedingly  complex.  The  first  method's  derivation  is  a 
reasonably  simple  task,  but  the  stability  of  the  resultant  filter  is  unpredictable. 
Also,  if  the  parameters  of  C  and  B  change  quickly,  the  linearization  will  be  invalid 
and  the  filter  will  diverge.  The  tuning  of  the  EKF  will  also  be  of  greater  complexity 
and  the  filter  will  be  computationally  more  burdensome  as  well.  Nevertheless,  if 
EKF  stability  for  this  model  could  be. shown  in  operation,  the  increased  tracking 
ability  of  this  filter  would  be  of  sufficient  worth  to  merit  its  use. 

One  method  of  incorporating  into  the  estimation  the  knowledge  that  the 
impulse  response  matrix  is  of  CB  form  without  linearized  or  nonlinear  filtering  is  to 
correct  the  Kalman  filter  or  WLS  estimate  so  that  it  is  of  the  same  form.  This 
process  is  referred  to  in  this  study  as  circulant  estimation  enhancement.  This  may 
be  done  by  a  variety  of  ad  hoc  methods.  The  method  used  here  "unwinds"  the 
estimate  matrix  through  permutations  of  each  column,  determines  the  B  diagonal 
elements  by  averaging  the  ratios  of  each  column  with  the  first  column  clement  by 
element,  and  then  averages  over  the  columns  to  determine  the  impulse  response. 
The  algorithm  is: 
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V  =[(1/Tn),(l/T2l),...(l/Tnl)] 

T1  =  [T,1.T,1....T111|T 

where  Tj  is  the  j’th  column  of  T  and  (Tj.)(j££  >s  the  k'th  column  of  the  corrected 
estimate.  It  is  important  to  emphasize  here  that  if  the  estimate  errors  of  T  are 
large,  the  above  procedure  will  result  in  an  even  poorer  estimate.  This  algorithm 
coupled  with  a  Kalman  filter  (which  would  use  in  the  next  recursion)  can 

therefore  produce  true  filter  instability.  Note  that  many  products  of  matrices  CB 
lie  near  T  (infinite  if  a  magnitude  constraint  is  not  placed  on  the  first  column  of  C 
or  the  trace  of  B).  An  Optimal  solution  for  C  and  B  (2n  parameters)  that  minimize 
a  measure  of  the  error  matrix  [CB  -  T]  could  be  computed  but  this  would  require 
solving  2n  nonlinear  equations  after  every  update.  (Another  CEE  algorithm  based 
on  unit  vectors  of  each  column  was  initially  tested  but  was  found  to  be  much  more 
sensitive  than  the  above  algorithm.) 
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CHAPTER  TWO:  CONTROLLERS 


2.1  Introduction 


With  the  plant  identified,  a  control  signal  can  be  fashioned  so  as  to  reduce 
the  vibration.  Typically,  optimal  control  theory  is  used  to  determine  a  control  law 
that  minimizes  a  performance  index.  The  control  law  is  just  the  functional 
relationship  between  the  control  that  minimizes  the  performance  index,  and  the 
plant  parameters,  filter  parameters,  plant  estimates,  and  measured  quantities  that 
comprise  the  information  state.  Usually,  for  a  discrete  plant  model  of  the  type 
employed  here,  the  performance  index  is  a  quadratic  of  the  output  vector  plus 
additional  terms  that  are  functions  of  the  input  vector  or  estimate  covariance. 
These  additional  terms  are  quadratics  attempting  to  discourage  certain  behaviors  or 
constraints  prohibiting  them.  The  performance  index  may  be  a  multistep  or  single 
step  formulation  of  all  these  terms  (output  vector  quadratic,  input  vector  functions, 
covariance  functions,  and  constraints). .  Depending  upon  whether  the  performance 
function  is  formulated  in  terms  of  the  expected  value  operator  or  the  identity 
operator,  is  single  step  or  multi-step,  and  on  what  terms  are  included,  certainty 
equivalent,  cautious,  dual,  or  sub— optimal  dual  control  laws  will  result  from  the 
minimization  of  the  performance  function. 

Three  additional  terms  are  added  to  the  vibration  quadratic  to  form  the 

performance  function  for  use  in  this  study.  They  are  a  control  quadratic  term, 

T  T 

(u  u)ip  a  control— rate  quadratic  term  (Au  Au)r9,  and  an  undesirable  control 

T 

mode  constraint,  u  'I' A.  Here,  u  is  the  mxl  control  vector.  An  is  the  control 
change  vector,  r^  and  r.,  are  weighting  factors,  U'  is  an  m xf  matrix  whose  columns 
represent  the  (  undesirable  control  modes,  and  A  is  a  (—  vector  of  undetermined 
Lagrangian  multipliers.  The  constraint  is  used  to  prevent  the  vibration  reduction 


control  from  initiating  rigid  body  motion.  For  this  study,  the  plant  model  has  been 
formulated  in  terms  of  a  quarter  revolution  per  cycle  in  vertical  vibration  only. 
Therefore,  the  only  rigid  body  mode  of  concern  is  vertical,  and  the  only  constraint  is 
that  the  control  has  zero  mean.  The  control  quadratic  term  is  included  to  encourage 
the  choice  of  low  magnitude  control,  since  the  accuracy  of  the  model  is  degraded  by 
large  blade  angles  of  attack.  (This  would  also  result  in  higher  blade  loads,  and 
possibly  actuation  difficulties).  The  control— rate  quadratic  is  included  so  as  to  limit 
actuator  power  requirements,  increase  the  accuracy  of  the  filter  model,  and  aid  in 
vibration  reduction  through  control  noise  filtering.  The  weighting  parameters  allow 
the  relative  importance  of  these  quadratic  terms  to  be  adjusted. 

One  interesting  modification  to  the  performance  function,  which  is  not 
included  in  this  research,  is  the  introduction  of  an  'internal'  control  rate  quadratic 
to  limit  the  control  variation  within  a  cycle.  This  would  be  useful  in  discouraging 
control  vectors  which  required  large  control  changes  within  a  cycle  and  therefore 
result  in  high  actuator  power  requirements.  This  can  be  formulated  by  including  a 
quadratic  of  the  differences  between  each  element  of  the  control  vector  and  its 
neighbors: 


(uT  G  u)  r.j  with  G  = 


2-1  0  •  •  •  0-1 

-1  2-1  0  •••  0 

0  -1  2  -1  0 

:  o  -l  2  ■  .  : 

0  '— 1 

-1  0  ••••  0  -1  •  2 


This  may  be  necessary  when  the  impulse  response  method  is  wind  tunnel  tested, 
especially  if  a  large  bandwidth  of  vibrations  is  to  be  reduced. 

In  the  investigation  for  this  project,  the  performance  function  was  written  as 
a  single  step  of  the  quadratic  terms  and  the  constraint.  Dual  control  laws  are  very 
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difficult  to  formulate  and  have  not  justified  their  complexity  sufficiently  in  previous 
studies  to  merit  inclusion  in  this  preliminary  investigation  of  impulse  response 
method  control.  They  do  have  many  interesting  properties,  however,  and  should  be 
investigated  in  further  studies. 

2.2  Cdoba.1  Certainty  Equivalent  Control  Law 
The  deterministic  performance  function, 


J  =  xTx  +  (uTu)i  J  +  (AuTAu)r,  +  uT*A  +  AT#Tu 


(2.1) 


where  all  quantities  without  subscripts  are  for  cycle  i+1,  yields  the  global  certainty 
equivalent  control  law  through  substitution  of  the  global  model  equation 


x  =  T  u  +  x 


o 


and  control  rate  vector  definition 


Au  =  u  —  U- 


(2.2) 

(2.3) 


where  is  the  cycle  i  control,  into  (2.1): 


rp  p  p  p  p  p  p 

J  =  xoxo  +  u  T  V+u  TiTu  +  xoTu  +  ui u(ri  +  r,>) 

rp  rp  rp  rp  T  T 

+  u  j  Uj  r0  -  u  uav,  -  u  ur9  +  u  M  + A  $  u 


Taking  the  first  variation  of  the  performance  function  yields 


(I)  =  duT  {TTx()  +  TTTu  +  u(r1  +  r.,)  -  u-  r9  +  4*  A} 
+  flAT{4fTu}  +  {uT*}  0A 
+  (uTT  1  +  x^  1  4-  u  (i‘j  +  r0)  —  u  J  r9  +  A^  4*} 


(2.4) 


(2.5) 


where  the  bracketed  terms  must  equal  zero  by  the  calculus  of  variations  when  <9J  is 
set  to  zero.  Therefore, 


[TTT  +  I(rj  +  r2)]  u  +  TTxq  -  m  r2  +  'M  =  0  (2.6) 

By  linear  algebraic  manipulation  the  control  law  is  found: 

Z  =  [TTT  +  I  (rL  +  i*2 )] 

u  =  —  Z  1  [TTxq  -  uj  r9  +  *  A] 

tfTu  =  o  =  -^TZ_1  [TTxq  -  Ujr9]  -  tfTZ_1\PA 

A  =  -[vpTz_1vp]-1  tfTZ-1  [TTxq  -  Ujr2] 

u  =  — Z— 1  [TTxq  -  Uj  r9]  +  Z"1  tf^Z-1*]-1  tfTZ-1  [TTxo  -  u.  r>2] 
u  =  {Z“1¥[¥TZ^1*r1¥T  -  1}  Z_1X  (2.7) 

where 

Z=[TTT  +  I(r1  +  r2)]  (2.S) 

X  =  [TTxq  -  Uj  r2]  (2.9) 

Because  of  the  similarities  between  the  performance  functions  used  in  this  study,  all 
the  optimal  control  laws  have  the  form  of  (2.7)  while  the  definitions  of  (2.8)  and 
(2.9)  for  the  matrix  Z  and  the  vector  X  are  different  in  each  case.  (This  form 
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results  from  the  orthogonality  constraint  term;  the  undesired  control  modes  are 


eliminated  through  the  use  of  a  Householder  transform.)  The  true  plant  values  xQ 
and  T  are  unknown,  therefore  the  filter  must  supply  estimates  to  be  used  in  (2.8) 
and  (2.9).  But,  since  the  control  for  cycle  i+1,  u,  is  computed  during  cycle  i  and  the 
output  y^  is  not  available  until  the  end  of  cycle  i,  new  a  posteriori  estimates  of  the 
parameters  are  unavailable  for  control  computation.  The  control,  therefore,  must 
use  the  cycle  i  a  priori  estimates,  xQ  and  Tj  in  (2.S)  and  (2.9).  Thus,  the  control 

derived 

Uj  j  =  {sr1*[*Tz-1*r1  y  - 1)  z-1x 

Z  =  [tT  T;  +  I  (rL  +  r,)]  (2.10) 

x  =  [TT  x  -  Uj  1‘2]  (2.11) 

is  optimal  with  respect  to  an  information  state  which  is  one  cycle  old.  The  resulting 
flow  of  measurements,-  estimates,  and  controls  between  the  plant,  estimator,  and 
controller  is  shown  in  Figure  2.1. 

2.3  Local  Certainty  Equivalent.  Control  Law 

By  substituting  in  the  local  model  relations 

Ax  =  T  Au  .  Ax  =  x  —  x. 

i 

a  (2  12) 

y  =  x  +  v  Au  =  u  —  u  \  > 

J  i 

for  x  and  u  the  performance  function  of  (2.1 )  becomes 
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rp  rp  T  T  T  T 

J  =  Ail  xT  Au  +  AuaTaXj  +  xjT  Au  +  xj  Xj 

4-  uT  Uj  +  Au^  Uj  +  uj  Au  i‘j  4-  (rj  4-  r.-,)  Au^Au  (2.13) 

4-  AuT  +  AT^TAu 

After  taking  the  first  variation  of  J  and  setting  it  equal  to  zero,  the  local  analogues 
to  (2.6)  are  obtained: 

[T  T  4-  I  (i'i  4-  iv,)]  Au  +  T  Xj  +  Uj  l'j  +  \&A  =  0 

(2.14) 

Au  =  0 


Here,  if  we  set 


Z  =  (TTT  +  I  (r,  +  r2)]  (2.15) 

X  =  (TT  Xj  +  u.  rj]  (2.16) 


algebraic  manipulation  yields  the  same  form  as  (2.7) 

Aui  +  1  =  {Z-1*  [tfT  Z_1  tf]-1  -  1}  Z-1  X 

Once  again,  to  obtain  the  quantities  to  use  in  (2.15)  and  (2.16)  we  must  use  the 
cycle  i  a  priori  estimate  IV  Furthermore,  since  the  best  estimate  of  Xj,  yj,  is 
unavailable  in  cycle  i  to  compute  the  cycle  i+1  control,  the  relation 


.  +  T.  Au. 

1  i  i 


x.  =  x- 

l  l 


(2.17) 
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must  be  used.  This  yields  the  local  certainty  equivalent  optimal  control  law: 


uj+1  =  Uj  +  {Z  1  #[*TZ  p  1  *T  -  1}  Z  1  X 


Z  =  [tT  T\  +  I  (i  j  +  r2)] 


X  =  [t7  +  T[  T,  Auj  +  Uj  rj 


(2.18) 

(2.19) 

(2.20) 


Note  that  this  law  feeds  the  measurement  noise  through  the  control  back  into  the 
plant.  This  can  significantly  affect  control  performance.  Since  the  control  feeds 
back  the  vibration,  this  is  a  closed-loop  control  law,  as  opposed  to  the  previous 
global  formulation  which  is  classified  as  an  open-loop  control  law. 


2.4  Global  Cautious  Control  Law 

The  global  cautious  control  law  is  obtained  via  the  minimization  of  the 
stochastic  performance  function 


J  =  E  {xrx  +  uTu(r1)  4-  AuTAu(r9)  +  uT®A  +  AT$Tu} 


(2.21) 


where  E  {  }  is  the  expected  value  operator.  By  inserting  in  (2.2)  and  (2.3)  this 
becomes 


T 


J  =  E  (xo‘xo  +  -TtTx„  +  “'T'T  u  +  xj  T  u) 

p  p  p  p  p  p  p 

+  u  u(r  1  +  r9)  +  u  .  u-r.,  -  u  miv,  -  u  j  ur.?  +  u  W  + A  f  u 


T  T 

i.  r 1  rr 


T  „ 


(2.22) 
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To  obtain  the  solution,  the  expected  value  term  must  be  expressed  in  terms  of  the 
available  information  state:  xQ  ,  Tp  P..  First,  consider  the  a  priori  covariance 

matrix  definition, 


Pj  =  E  {(71  -  tfi)T  (^-dj)} 


II 

■3a 

X 

T. 

1 

L  °i 

l 

The  covariance  matrix  can  be  partitioned  as 


p12j 

K 

p 

r22j 

Pll;  =  E  {(X0:  “XOj)T  (X0;  -XOj)i  =  E  ^ 


P19.  =  E  {(x0  -  x  )T  (Tj  -  T  )}  ee  E  { (  T  <T  } 
“1  ii  '  i  i 


P-n.  =  E  {(T  - T;)T  (X  -x  )}  =  E  { CTT  Cx  ) 
I  i  i  i  '  i 


P.»  =  E  {Cl-,  -  T;)T  (T, -T  )}  =  E  {  CtT  <t  } 
i  i  i 

C  =  x  —  x 
xi  °i  °i 


CTj  =  Ti  -  Ti 


(2.23) 


(2.24) 


(2.25) 
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and  the  cycle  i+1  and  i  plants  are  related  by  the  filter's  random  walk  model 


+  w 


x- 


T  =  Tj  +  wT 


(2.26) 


with  E  {w^  wx  }  =  Wn,  E  (w^  wT  }  =  0^,  E  {w!£  wx  }  =  On,  and 
T 

E  {w^  (On  is  a  null  n— vector.  Note  that  the  process  noise  here  is 

i  i 

stationary.).  Using  (2.25)  and  (2.26)  the  cycle  i+1  quantities  x  ,  T  can  be  related 
to  the  cycle  i  information  state, 


*„  =  VVw*i 


T  —  I'i  Cm  +  Wrp 
i  i 


(2.27) 


and  the  expected  value  term  of  (2.22)  can  be  solved: 


P  f  T  ,  r  r-T  -  -T  ,  ,  -T 

L  x  =  L  x„  x„  -x„  (  +  w 

1  O  0J  1  0  o.  O-^X.  O-  X- 

II  11  11 

+  CT  C  -  CT  X  -  CT  W  (2.28) 

xi  xi  xi  °i  xi  xi 

,  T  T  >  .  T  —  i 

+  w  w  —  w  C+w  x 

X-  X-  x.  \x.  X-  O  j 

11  11  II 

Since  the  estimate  is  a  constant,  it  passes  through  the  expected  value  operator. 
Therefore,  with  xQ  an  unbiased  estimate  and  a  zero  mean  random  vector. 


E  el  U  =  ;l K  IS.)  = 0 

ii  i  i 
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E  wx  }  =  xo  E  {wx  }  =  0 
i  "  i  i  '  i 


and  the  random  walk  statistically  independent  of  the  estimate  error, 


E  {wT  C  }  =  0 

"i  '  i 


(2.28)  simplifies  to 


E  {xT  x  }  =  xT  x  +  E  (CT  C  }  +  E  fwT  w  } 

1  O  0;  OO-  lsXj  \\jJ  1  Xj  XjJ 


E  lxo  xo)  =  xoixoj  +  piij +  wn 


(2.29) 


Similarly. 


I?  I  TrpT  1 

E  {u  r  xo} 


=  uT  E  {TTx0( 

=  uT  E  {(T,  -  CT  +  )T  (x  -  C  +  w  )} 

i  i  i  '  i  J  i 

=  uT  (TT  x  -T1  C  +  tT  w 

v  1  0-  1  SX-  1  X- 

1  1  1 

,T  *  ,T  -  ,T 

+  ^T.  <x.  ^T.  xo.  ^T.  wx. 
ii  it  ii 

T  T  T  — 

+  wT  w  -  wT  C  +  wT  x  } 

i  '  i  i  '  i  i  i 

=  uT  (T[  xq  +  E  Cx  }  +  E  {wj  w  }) 
i  i  '  i  i  '  i 


r  /  T,,,T  , 

E  {u  1  xo} 


ljT  (  l  |  Xq  +  1*2! 


(2.30) 


and  by  the  same  methodology, 
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E  {  uTTTT  u}  =  uT(t]  Tj  +  P99  +  W99)  u 

"“i 


(2.31) 


Equations  (2.29),  (2.30),  (2.31),  and  (2.30)  transpose  define  all  the  terms  in  (2.22); 
therefore, 


J  =  <*0.  *0.  +  Pll.  +  Wll>  +  “T  <TI  *o.  +  p21.> 

111  11 

+  Ti  +  Fl->  )  “  +  “T  (T?  T,  +  p»  +  W.„)  u  (2.32) 

i  “i  ““i 

/T*  T1  rp  rp  rp  rp  rp 

+  u  u(i'j  +  r9)  +  u  r  Uji9  -  u1  U|i9  -  u|ur9  +  u 1  'I'A  +  A 1 


Taking  the  first  variation  yields 


dJ  = 


°"T  [<TI  *0j  +  P21j>  +  <tT  Ti  +  1*22.  +  W22>  u 

+  u(rj  +  r9)  -  Uj(r2)  +  tfA] 

+  dA^  [^u]  +  [uT*]  dA 

+  i<xo  h + Fi-> )  + uT  *T!  h + +  w,.,) 

i  “i  ““i 

+  uT(ri  +  r9)  -  u|(r9)  +  A^'V*']  cTu 


Setting  this  equal  to  zero  results  in  the  two  equations: 


I’l'i  Tj  +  I’.,,  +  W.„  +  Id  .  +  r.,)]  u  =  ['1'T  x  +  I*  -  u-( r., )|  +  *A 
““i  “  i  “  i 

(2.33) 


* 1  u  =  0 


> 


These  are  the  cautious  analogues  to  the  (2.6)  equations  and  have  the  solution: 


Uj  .  ,  =  (Z  1  H  (*T  Z  1  *]  1  *T  -  1}  Z  1  X 

Z  =  [Tj  Tj  +  I(r,+r,)  +  F.„  +  W,,]  (2.34) 

i 

X=(T'[;0“ui(r2)  +  P2lJ  (2-33> 

This  is  the  global  cautious  control  law.  Note  the  difference  between  these  equations 
and  (2.10)  and  (2.11).  The  addition  of  the  covariance  and  process  noise  term  in  the 
Z  matrix  formulation  acts  like  increased  control  quadratic  weighting.  This 
discourages  large  control  vector  magnitudes  when  the  Filter  is  uncertain  about  the 
impulse  response  matrix  estimate.  The  cross  covariance  vector  P.-^  in  (2.35)  acts 
the  opposite  of  control  rate  quadratic  weighting  encouraging  change  most  heavily  in 
those  elements  of  the  control  vector  that  are  associated  with  large  uncertainties  in 

tTV 

2.5  bocal  Cautious  Control  Law 

By  taking  the  first  variation  of  the  stochastic  performance  function  of  (2.21) 
with  the  local  model  substitutions  (2.12)  and  setting  it  equal  to  zero,  the  following 
two  equations  are  obtained: 


43 


As  was  found  for  the  global  cautious  control  law, 

E  {TT  T}  =  T[  Tj  +  P99  +  W99 

(where  F99  is  now  the  entire  covariance  matrix  since  the  local  model  Kalman  filter 
““i 

only  estimates  T-  and  not  x  )  whereas  the  expected  value  term  on  the  right  side  of 

1 

(2.36)  has  not  been  encountered  before.  Immediately  it  comes  to  mind  to  replace  x. 

in  this  term  with  y-  since 
J  i 

E  {*()  =  », 

because  the  measurement  noise  is  zero  mean.  However,  the  measurement  vector  y. 
is  not  available  to  the  end  of  the  cycle  i,  while  the  control  for  cycle  i+1  must  be 
computed  during  cycle  i.  Therefore,  the  quantity  x^  must  be  related  to  the  previous 
measurement: 

x  =  x-  ,  +  T-  Au- 
li—l  ii 

x.  =  y.  ,  —  v-  ,  +  T-  Au- 
i  Ji— 1  l—l  i  l 

E  {Tr  x.}  =  E  {TT  (y-  .  -  v.  .  +  T.  Au-} 

i  iJ  1  VJi— 1  i—l  i  iJ 

from  (2.27) 
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E  {TT  Xj}  =  E  {(1'T  -A  +  )(y.  -  v.  ,  +  (T.  -  CT  )  Au.)} 

•  i  i  i 

Assuming  the  estimate  errors  are  uncorrelated  with  the  measurement  noise 

E  (T  xj  =  TT  yj_1  +  (TT  Tj  +  P,2  )  AUj 

•  and  the  solution  to  (2.36)  is 

ui+i =  ui +  <z_1  *  i*"1  z_1  *r‘ 4,1  - 1)  z'1  x 

Z  =  [xT  T.  +  I  (r,+r2)  +  P22  +  W,,]  (2.37) 

•  X  =  [TT  yw  +  T{  T;  Auj  +  Uj(r,)  +  P22  AUj) 

This  control  law  contains  covariance  and  process  noise  matrices  not  present  in  the 
local  certainty  equivalent  control  law  in  (2.20).  The  role  of  the  covariance  matrix  in 
this  equation  is  similar  to  that  of  the  control  quadratic  weighting.  Through  the  Z 
matrix  equation,  the  uncertainty  acts  to  decrease  the  change  in  the  control,  while 
through  the  X  vector  equation  the  uncertainty  acts  to  excite  changes  in  those  modes 
of  the  control  vector  whose  previous  changes  lied  close  to  the  larger  principal  axes  of 
the  covariance  ellipsoid. 

2.6  Performance  of  Ccrtnintv  Eouivnlcnt  and  Cautious  Controls 


As  pointed  out  by  Bar— Shalom  [50].  cautious  and  certainty  equivalent 
controls  are  best  suited  for  different  control  problems.  Cautious  policies  seek  to 
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avoid  the  cost  penalties  associated  with  control  of  an  unknown  system.  Thus,  they 
function  best  on  problems  where  the  initial  uncertainty  is  high  and  control  is 
unlikely  to  result  in  better  future  estimates.  Certainty  equivalent  controllers,  while 
they  do  not  actively  probe  the  system  like  dual  controllers,  are  more  vigorous  and 
thus  incur  both  higher  uncertainty  costs  and  better  estimates.  If  the  system  varies 
slowly  enough  and  useful  information  can  be  extracted  above  the  noise  of  signal, 
these  controllers  will  improve  their  estimates  and,  therefore,  their  control.  If  the 
initial  estimates  are  poor  and  the  time  horizon  is  sufficiently  short,  the  exacerbation 
of  the  cost  through  the  poor  control  will  not  be  compensated  by  the  later 
improvements  in  control  due  to  better  identification. 

2.7  Inversion  of  the  Z  Matrix 

One  significant  problem  with  the  optimal  control  laws  previously  formulated 
is  the  inversion  of  the  Z  matrix.  This  matrix  is  mxm  (where  m  is  the  size  of  the 
control  vector,  which  is  typically  equal  to  n,  the  size  of  the  measurement  vector) 
and  thus  requires  m^  operations.  For  a  quarter  revolution  cycle  and  vibrations 
harmonics  up  to  1614  controlled,  m  must  be  at  least  an  Sxl  vector  by  the  Nyquist 
requirement.  More  realistically,  the  system  should  have  at  least  four  times  this 
sampling  frequency  to  yield  proper  control  (m  =  32).  A  matrix  of  this  size  could 
not  be  inverted  in  real  time  with  present  methods.  This  problem  may  be  solved  in 
several  ways.  With  the  advent  of  parallel  computing  and  the  increasing  speed  of 
processors,  in  the  near  future,  matrices  of  this  size  may  be  invertible  in  real  time. 

Iterative  routines  exist,  such  as  the  Pan-Rief  algorithm  [60],  which  can 
quickly  find  a  matrix  inverse  if  supplied  with  a  close  approximation  (such  as  the 
previous  cycles  Z  ^).  The  controller  algorithm  can  be  adjusted  (for  the  cautious 
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case  this  would  require  re— derivation)  so  that  the  lag  time  between  the  obtaining  of 
Tj  and  its  use  in  the  control  is  greater  than  the  present  one  cycle.  (The  estimate  xQ 

could  still  be  used  in  computing  the  cycle  i+1  control.)  This  method  would  lessen 
the  vibration  reduction  obtained  by  the  regulator  in  proportion  to  the  quickness 
with  which  T  changes. 

Another  method,  which,  while  primitive,  may  work  sufficiently  well  in  the 
field  would  have  the  local  Kalman  filter  identify  T  *  (provided  m=n).  Then,  if  a 
certainty  equivalent  policy  is  used,  the  constraint  condition  is  ignored,  and  no 
control  or  control— rate  weight  is  used  (which  as  shall  be  seen  with  the  local  model 
produces  good  control),  then  the  control  law  becomes 


u  =  „i  —  ( T  Vi-, 

where  (T  ^).  is  the  a  priori  estimate  of  impulse  response  matrix  inverse.  This  is  a 
form  of  adaptive  inverse  control  [61].  Note  that  no  matrix  inversion  is  required. 

If  the  zero  mean  control  constraint  is  dropped,  the  optimal  control  law  (for 
both  certainty  equivalent  and  cautious  policies)  is 
u  =  -  Z_1  X 

which  is  the  solution  of 

Z  u  =  -X 

The  optimal  control  u  can  be  found  in  n'^/3  operations  if  it  is  not  necessary  to 
obtain  Z  Thus,  a  66%  reduction  in  the  number  of  operations  can  be  achieved.  If 


the  performance  function  is  augmented  with  an  orthogonality  quadratic 


and  the  weight  1^  is  large,  then  the  resulting  control  will  approach  the  constrained 
control  law.  The  Z  matrix  will  then  contain  an  orthogonality  term: 

Z  =  (T?  Tj  +  I  (r,+r2)  +  *4’T(r4)] 

As  discussed  in  the  previous  chapter,  the  impulse  response  matrix  is  the 
product  of  a  circulant  matrix  and  a  diagonal  matrix.  If  the  estimate  could  be 
decomposed  into  these  two  parts,  then  the  inversion  of  the  nxn  certainty  equivalent 
Z  matrix  can  be  accomplished  in  0  (n  log.-,  n)  operations  if  the  control  and 
control— rate  weights  are  zero. 

z  =  [tT  T] 

Tj  =  Cj  Bj  ,  Cj  =  nxn  circulant  matrix 
Bj  =  nxn  diagonal  matrix 

Z  =  B  CT  C.  B. 

■  1111 

Since  B  is  diagonal  it  is  invertible  in  n  operations.  Also,  because  the  multiplication 
of  two  circulants  yields  a  circulant  and  any  circulant  can  be  inverted  in  0  (n  log.,  n) 
operations  via  FFT  techniques,  the  computation  of 

zr1  -  B  _l  (CT  c.)_1  B._1 

can  be  completed  in  O  (n  log.,  n).  For  n=T2  this  can  result  in  99%  fewer  operations 
than  if'*  full  matrix  inversion.  Note  that  the  zero  mean  control  constraint  can  still 


lie  used  (as  opposed  to  the  previous  method). 


2.8  Output  Constraint 

T 

A  control  constraint,  u  ^  =  0  ,  was  used  in  the  previously  defined 
performance  functions  to  prevent  control  modes  which  would  initiate  rigid  body 
motion.  Since  the  azimuthal  blade  effectiveness  is  dependent  on  the  flight 
condition,  the  elements  of  will  vary  during  flight.  Thus,  any  fixed  control 
constraint,  like  that  used  in  this  paper,  will  be  only  approximately  true  at  any  time 
during  flight.  Another  method  to  prevent  control  induced  rigid  body  motion 
introduces  an  output  constraint, 

T 

x'l'  =  0 

for  the  certainty  equivalent  performance  function,  or 

E(xTtf)  =  0 

for  the  cautious  performance  function.  After  minimization,  the  control  laws 

ui+1  =  [Z”1tT$  (^TjZ^T^r1  *T  Tj  -  I]  Z-1  X  (global) 

uj+1  =  Uj  +  [Z_1tT«  (^TjZ^tT*)-1  Tj  -  I]  Z_1  X  (local) 


result.  These  are  the  equivalents  of  (2.7)  and  (2.18).  The  terms  Z  and  X  are  the 
same  as  those  derived  in  Sections  2.2,  2.3,  2.4,  and  2.5  depending  on  the  type  of 


control  law. 

The  prime  disadvantage  of  the  output  constraint  is  that  the  efficacy  of  the 
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control  (and  the  freedom  from  induced  rigid  body  motion)  is  dependent  on  the 
accuracy  of  the  impulse  response  matrix  estimate.  For  this  reason,  it  was  rejected 
for  this  investigation. 


cycle  i-1 


cvcle  i 


cycle  i+1 


Figure  2.1  Time  sequence  for  the  plant,  estimator,  and  controller 
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CHAPTER  THREE:  RESULTS.  CONCLUSIONS.  AND  RECOMMENDATIONS 

3.1  Introduction 

The  estimation  and  control  algorithms  developed  in  the  previous  two 
chapters  were  tested  via  simulation  with  the  computer  model  described  in  Section 
1.5.  Several  tests  were  designed  to  evaluate  the  performance  of  the  filters  and 
controllers.  Open  loop  esitmation  simulations  were  conducted  first  to  establish  an 
understanding  of  filter  behavior.  The  closed  loop  regulators  were  then  tested  and 
analyzed.  These  results  will  be  discussed  but  not  presented  here  —  for  more  detail 
please  see  Reference  63. 

3.2  Simulation  Results 

The  Kalman  filter  regulators  consistently  yielded  between  89  and  93% 
vibration  reduction  in  the  simulations  executed.  Due  to  their  superior  noise 
suppression  and  better  handling  of  uncontrolled  vibration  changes,  the  global 
Kalman  filter  regulators  performed  slightly  better  than  their  local  counterparts. 

The  invariable  global  and  local  regulators  produced  fairly  good  results  on 
most  of  the  tests;  in  situations  where  there  is  little  change  in  the  plant,  these  can 
outperform  the  adaptive  regulators.  Since,  for  global  Kalman  filters,  tunings  that 
accentuate  impulse  response  matrix  tracking  diminish  uncontrolled  vibration 
tracking  and  vice  versa,  it  is  highly  desirable  to  base  the  tuning  and  filter  model 
used  on  the  types  of  plant  changes  most  likely  to  be  encountered.  If  the  plant 
changes  are  concentrated  in  either  one  of  T  or  xQ,  the  optimal  tuning  regulator  will 
in  essence  be  invariable  to  the  other.  This  fact  may  be  used  to  reduce  the 
computational  burden  of  the  filter  and  controller  without  degrading  regulator 
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performance.  For  example,  if  the  impulse  response  matrix  changes  very  little  under 
typical  flight  conditions,  then  a  T  invariable  global  regulator  could  be  used 
providing  quick  response  to  changes  in  uncontrolled  vibration  with  a  simpler  filter 
and  an  optimal  control  solution  not  requiring  inversion  of  Z. 

It  is  interesting  to  note  that  local  regulators  performed  well  when  confronted 
with  quick  changes  in  uncontrolled  vibration  even  though  this  change  violates  their 
formulations.  This  is  at  first  confusing  since  open  loop  local  estimators  were  weak 
in  this  area.  But,  the  control  variations  produced  by  the  changes  in  the  plant 
resulted  in  increased  identification  prowess  allowing  for  quick  recovery.  When  the 
plant  is  changing  slowly,  the  Kalman  filter  was  found  to  track  successfully.  Thus, 
even  with  small  control  variations,  useful  information  can  be  obtained  about  a  noisy 
slowly  drifting  plant,  permitting  continuous  vibration  reduction.  The  closed  loop 
regulator  therefore  appears  to  have  avoided  in  practice  its  paradoxical  kernel:  for  a 
stationary  plant,  full  rank  control  is  necessary  for  complete  identification  yet 
complete  identification  implies  rank  one  control. 

The  batch  weighted  least  squares  filter  regulators  typically  reduced 
vibrations  S3  to  85%  when  optimally  configured.  The  local  WLS  regulators 
performed  better  than  its  global  counterpart  due  to  the  level  of  inaccuracy  in  the 
uncontrolled  vibration  estimate.  For  all  but  the  random  walk  test,  the  optimal 
batch  size  was  twelve.  (In  the  random  walk  test  the  plant  does  not  move  as  quickly 
away  from  its  current  state  as  in  the  other  tests;  therefore,  a  larger  batch  size 
improves  identification.)  This  size  exemplifies  the  classic  principle  of  filter 
sensitivity:  high  sensitivity  causes  noise  to  be  mistaken  for  information;  low 
sensitivity  causes  information  to  be  mistaken  for  noise.  Neither  in  open  loop  or 
closed  loop  does  the  WLS  filter  identify  as  accurately  as  the  Kalman  filter.  In  the 


open  loop  simulations  this  is  largely  due  to  measurement  noise.  The  control 
constraint  compounds  this  problem  in  the  closed  loop  simulations. 

Certainty  equivalent  control  performance  was  degraded  by  the  inclusion  of 
any  control  weight.  The  use  of  a  small  control— rate  weight  resulted  in  greater 
vibration  reduction  for  global  regulators  by  limiting  the  control  response  to  estimate 
errors  (thus  acting  to  filter  out  measurement  noise).  For  local  regulators,  however, 
the  control  must  remain  free  to  follow  changes  in  uncontrolled  vibration  not 
represented  in  the  plant  model  quickly. 

Cautious  regulators  even  after  retuning  do  not  perform  as  well  as  certainty 
equivalent  regulators.  While  they  approach  these  levels  of  vibration  reduction,  they 
are  still  sluggish.  Caution  does  yield  some  smoothing  in  behavior,  but  this  can  be 
achieved  through  control-rate  weight  and  thereby  the  increased  complexity  of  the 
cautious  algorithms  can  be  avoided. 

Two  methods  were  explored  to  prevent  the  control  constraint  induced 
covariance  instability  in  the  Kalman  filter  regulators.  The  first,  random  probing 
insertion,  violates  the  control  constraint  deterministically  but  not  statistically.  This 
caps  the  covariance  growth  but  degrades  the  vibration  reduction  performance. 
Philosophically,  probing  is  poorly  grounded  since  it  cures  a  filter  ailment  through 
the  corruption  of  the  optimal  control.  The  second  method  advanced,  covariance 
trace  bounding,  alters  the  filter  algorithm  directly  to  prevent  covariance  instability. 
This  method  is  justified  by  a  simple  truism:  the  difference  between  a  filter  estimate 
and  the  true  value  will  always  be  well  bounded.  Covariance  trace  bounding  was 
found  to  not  degrade  vibration  reduction  performance  (unless,  of  course,  the  bound 
was  unrealistically  low).  When  this  method  was  used  with  high  process  noise 
variance! s)  the  filter  attained  nearly  fixed  covariance  eigenvalues.  It  may  be 
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possible  that  the  filter  algorithm  could  use  a  fixed  covariance  matrix  simplifying  the 
computations  of  the  Kalman  gain  vector  and  the  estimate  updates. 

Any  increase  in  measurement  noise  reduces  identification  accuracy  and 
thereby  vibration  suppression.  For  the  Kalman  filter  regulators,  the  degradation  in 
vibration  reduction  is  approximately  linear  in  the  measurement  noise  standard 
deviation  over  the  range  tested.  When  the  noise  had  standard  deviation  0.1  (ten 
times  the  baseline)  the  global  Kalman  filter  regulator  produced  the  best 
performance  with  83%  reduction  while  the  local  WLS  filter  regulator  had  the  worst 
with  only  68%  vibration  reduction.  At  this  level  of  noise,  the  local  filter  regulators 
perform  poorly  since  they  feed  the  noise  directly  into  the  control. 

The  circulant  estimation  enhancement  algorithm  failed  to  improve  the 
vibration  reduction  performances  when  added  onto  the  filter  algorithms  even  though 
it  appeared  to  benefit  estimation  performance  in  the  open  loop  simulations.  This 
most  likely  resulted  from  losses  in  accuracy  along  some  axes  of  the  error  ellipsoid 
while  gains  occurred  along  others.  In  some  cases  if  the  filter  estimates  were  poor, 
the  CEE  algorithm  would  exacerbate  the  estimate  errors.  This  resulted  for  a 
number  of  open  loop  and  closed  loop  Kalman  filter  simulations  in  true  instability. 
(Since  the  WLS  algorithm  is  non— recursive,  it  is  always  stable.)  This  instability 
proved  to  be  removable  through  re— tuning.  The  CEE  global  Kalman  filter  regulator 
did  significantly  outperform  all  other  regulators  in  high  noise  testing.  When  a 
measurement  noise  standard  deviation  of  0.1  was  used  an  additional  3%  of  vibration 
reduction  was  obtained.  While  the  CEE  algorithm  presented  in  this  report  is  clearly 
unacceptable  due  to  its  instability,  many  other  ad  hoc  (possibly  sub-optimal) 
algorithms  could  be  formulated.  One  of  these  may  yield  better  tracking  as  well  as 
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noise  suppression  without  either  instability  or  undesirable  complexity.  An  extended 
Kalman  or  full  nonlinear  filter  could  also  be  employed. 

For  greater  detail  on  the  open  and  closed  loop  results  please  see  Reference  63. 

3.3  Summary.  The  Impulse  Response  Method 

A  new  method  of  tackling  the  helicopter  vibration  problem  was 
demonstrated.  The  plant  is  characterized  in  terms  of  a  discrete  sequence  of 
uncontrolled  output  and  a  matrix  of  impulse  responses  in  order  to  fashion  a  discrete 
control  sequence  to  minimize  a  performance  function  quadratic  in  the  output.  If  the 
uncontrolled  output  is  predictable  and  the  plant  dynamics,  and  therefore  impulse 
response,  are  nearly  stationary,  then  the  output  can  be  controlled  with  the 
assistance  of  an  estimator.  The  periodicity  of  the  helicopter  plant  yields  the 
predictability  of  the  uncontrolled  output  and  sets  the  cycle  length  of  the  regulator. 
For  the  control  of  vertical  vibration  the  cycle  need  only  be  one  quarter  of  a  blade 
revolution.  If  a  nearly  constant  control  condition  is  applied,  the  impulse  responses 
in  the  matrix  "wrap-around"  producing  a  full  matrix.  This  assumption,  which  is 
similar  to  the  quasi-static  assumption  of  higher  harmonic  control,  produces  a 
circulant  matrix  for  systems  with  constant  dynamics  and  a  circulant— diagonal- 
product  matrix  for  systems  with  periodic  input  coefficients.  Circulant  matrices  can 
be  diagonalized  via  Fast  Fourier  and  inverse  Fast  Fourier  transforms;  therefore,  the 
impulse  response  matrix  formulation  is  directly  related  to  both  discrete  time  and 
frequency  (HHC)  domains.  Both  relationships,  and  the  impulse  response  method 
itself,  merit  a  theoretical  investigation  in  their  own  right. 
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3.4  .Summary:  Kalman  and  Weighted  Least  Squares  Filtering 

Global  and  local  Kalman  filters  were  derived  to  estimate  the  plant 
.  parameters.  The  resulting  algorithm  was  computationally  simplified  by  the 
assumption  of  equivalent  process  and  measurement  noise  covariances  for  each  row  of 
the  plant  variable  matrix.  Global  and  local  batch  weighted  least  squares  (WLS) 
algorithms  were  also  found.  In  regulator  application  the  control  constraint 
confounded  the  operation  of  the  algorithms  by  preventing  the  matrix  to  be  inverted 
from  attaining  full  rank.  This  was  solved  by  adding  pseudo— random  noise  to  the 
control  vector  so  that  the  control  constraint  was  kept  in  a  statistical  sense.  A  more 
elegant  means  to  attain  full  rank  might  be  to  impose  an  additional  constraint  on  the 
plant  variable  matrix. 

It  was  demonstrated  theoretically  and  experimentally  that  the  stability  of 
the  Kalman  filter  covariance  equation  required  a  matrix  of  controls  to  have  full 
rank.  Furthermore,  it  was  shown  that  the  number  of  bounded  covariance 
eigenvalues  was  equal  to  the  rank.  This  simply  means  that  if  information  is  not 
received  by  the  filter  along  an  axis  of  the  covariance  ellipsoid  then  the  uncertainty 
along  this  axis  will  grow  due  to  the  process  noise  covariance.  This  presented  a 
problem  to  closed  loop  Kalman  filter  performance:  the  use  of  the  control  constraint 
guarantees  that  the  control  matrix  will  be  less  than  full  rank.  Three  solutions, 
covariance  trace  bounding,  insertion  of  random  noise,  and  filter  turn— off/on,  were 
offered.  The  first  two  of  these  were  incorporated  into  the  filter  subroutines.  (It  is 
difficult  to  formulate  a  condition  to  allow  the  filter  to  turn— on  again;  also,  this 
solution  would  disregard  information  along  the  axes  that  the  control  is  actuating.) 

The  cireulant/diagonal  structure  advanced  for  the  impulse  response  matrix 
theoretically  permits  identification  in  two  cycles.  The  filtering  process,  however, 
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must  be  considerably  more  complex  since  the  model  is  nonlinear  in  these 
parameters.  A  regulator  based  upon  this  model  would  track  changes  in  the  plant 
with  extreme  agility  although  its  stability  would  be  questionable.  The  simplest 
identification  algorithm  for  this  nonlinear  problem  would  be  an  extended  Kalman 
filter. 

In  this  research  the  circulant/diagonal  structure  of  the  impulse  response 
matrix  could  be  utilized  without  performing  full  nonlinear  filtering  by  a  method 
called  circulant  estimation  enhancement  (CEE).  Here,  the  structure  was  imposed 
on  the  estimate  supplied  by  the  filter  through  an  ad  hoc  algorithm.  Many  CEE 
algorithms  could  be  derived  including  an  optimal  one  which  would  require  solving 
2n  nonlinear  equations  after  each  update. 

3.5  Summary:  Control  Laws 

Four  optimal  control  laws  employing  the  filter  estimates  were  derived:  the 
global  certainty  equivalent,  the  local  certainty  equivalent,  the  global  cautious,  and 
the  local  cautious.  These  were  optimal  with  respect  to  an  information  state  which 
was  one  cycle  old.  All  the  performance  functions  contained  a  zero  mean  control 
constraint  to  prevent  vertical  rigid  body  motion.  The  four  control  laws  derived 
have  an  identical  form  because  the  constraint  condition  results  in  the  removal  of  the 
undesired  modes  via  a  Householder  transform.  This  constraint  was  an 
approximation  since  the  actual  lift  generated  by  the  four  blades  is  a  function  of 
blade  azimuth.  However,  this  approximation  allowed  the  constraint  to  be 
independent  of  the  estimate  accuracy  and  was  therefore  considered  desirable. 
Considering  the  accuracy  of  the  filter  estimates  maintained  during  the  simulations 
presented,  an  output,  constraint  dependent  on  the  estimates  may  present  no 


problems  to  regulator  performance  whatsoever.  In  addition  to  restricting  rigid  body 
motion  more  accurately,  the  output  constraint  would  not  result  in  rank  deficient 
control  due  to  the  -andom  fluctuations  in  the  estimates.  Thus,  a  'natural'  probing 
would  be  present  and  the  covariance  equation  of  the  Kalman  filter  would  be 
stabilized.  If  WLS  filtering  was  used,  the  control  batch  matrix  would  be  of  full  rank 
despite  the  constraint. 

The  optimal  control  laws  presented  require  the  inversion  of  the  mxm  matrix 
Z.  For  multi— axis  or  high  bandwidth  vibration  control  this  could  not  be  performed 
in  real  time  onboard  as  it  requires  0  (n  )  operations.  Several  methods  to  overcome 
this  problem  were  advanced  including  parallel  computing  algorithms,  iterative 
routines,  adaptive  inverse  control,  and  employment  of  the  circulant/diagonal  form. 
In  the  last  of  these  methods,  the  inversion  is  possible  in  0  (nlog9n)  operations. 

3.6  Recommendations 

At  this  preliminary  stage  of  research  into  the  application  of  the  impulse 
response  method  to  the  helicopter  vibration  problem,  many  different  avenues  of 
investigation  are  open  and  worthy  of  further  study.  Those  of  greatest  import  to  the 
furtherance  of  the  method  are: 

(1.)  Continued  Theoretical  Development  of  the  Impulse  Response  Method: 
The  method,  independent  of  application,  needs  to  be  firmly  and 
completely  connected  to  the  body  of  control  theory.  Circulant 
relationships  should  be  among  the  tools  exploited  in  this  process. 
Application  to  different  types  of  systems  (constant  coefficients, 
periodic  coefficient  )  should  be  examined. 


(2.)  Demonstration  of  the  Method:  The  method  should  be  tested  in  the 
control  of  a  one  degree  of  freedom  dynamical  system  with  a  periodic 
forcing  function  and  periodic  coefficients.  This  should  be  compared  to 
results  of  a  higher  harmonic  control  approach.  An  extended  Kalman 
filter  impulse  response  regulator  using  a  circulant  structure  should 
also  be  tested. 

(3.)  Further  Development,  of  Helicopter  Multi— Axes  Models:  These 
models  should  be  developed  in  the  complete  context  of  blade  dynamics 
and  the  rotor  flow  field  so  as  to  be  understandable  by  the  general 
helicopter  community. 

These  recommendations  are  designed  to  remedy  the  skepticism  and  confusion  of  the 
helicopter  community  over  this  new  control  approach.  Other  areas  of  investigation 
such  as  the  matrix  inversion  problem,  output  constraints  in  the  performance 
function,  and  plant  models  that  are  not  based  on  the  nearly  constant  control 
assumption  are  certainly  also  of  importance.  But,  advances  in  these  areas  will  not 
appear  to  be  of  sufficient  worth  to  pique  the  insular  interests  of  the  helicopter 
community.  It  is  hoped  that  by  expanding  its  theoretical  base,  demonstrating  its 
application,  and  suggesting  algorithms  for  actual  helicopter  implementation,  the 
impulse  response  method  will  be  received  with  increased  interest. 
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